Testing my shiny new cosmological emulator.


In [1]:
from pearce.emulator import OriginalRecipe, ExtraCrispy
from pearce.mocks import cat_dict
import numpy as np
from os import path

In [2]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()

In [3]:
training_data = '/u/ki/swmclau2/des/PearceTrainerTest1.hdf5'

em_method = 'gp'
split_method = 'random'

In [4]:
a = 1.0
z = 1.0/a - 1.0
print z


0.0

In [5]:
fixed_params = {'HOD':4, 'r':1.07508818}#, 'r':0.18477483}

In [16]:
n_leaves, n_overlap = 2, 1
emu = ExtraCrispy(training_data, n_leaves, n_overlap, split_method, method = em_method, fixed_params=fixed_params)


['ombh2' 'omch2' 'w0' 'ns' 'ln10As' 'H0' 'Neff']
[[  2.26832500e-02   1.14059800e-01  -8.16597200e-01   9.75589000e-01
    3.09291800e+00   6.33656900e+01   2.91875000e+00]
 [  2.24507100e-02   1.17257600e-01  -1.13444000e+00   9.76522700e-01
    3.14987200e+00   7.30997200e+01   3.17375000e+00]]
emu = OriginalRecipe(training_data, method = em_method, fixed_params = fixed_params)

In [17]:
emu._ordered_params


Out[17]:
OrderedDict([('ombh2', (0.022450709999999999, 0.022683249999999999)),
             ('omch2', (0.11405979999999999, 0.1172576)),
             ('w0', (-1.1344399999999999, -0.81659720000000002)),
             ('ns', (0.97558899999999993, 0.97652269999999985)),
             ('ln10As', (3.0929180000000001, 3.1498720000000002)),
             ('H0', (63.365690000000001, 73.099719999999991)),
             ('Neff', (2.9187500000000002, 3.1737500000000001)),
             ('z', (0.0, 0.0))])

In [8]:
print emu.y


[ 0.08626049 -0.08626049]

In [9]:
emu.x


Out[9]:
array([[  2.26832500e-02,   1.14059800e-01,  -8.16597200e-01,
          9.75589000e-01,   3.09291800e+00,   6.33656900e+01,
          2.91875000e+00,   0.00000000e+00],
       [  2.24507100e-02,   1.17257600e-01,  -1.13444000e+00,
          9.76522700e-01,   3.14987200e+00,   7.30997200e+01,
          3.17375000e+00,   0.00000000e+00]])

In [10]:
emu._ordered_params.keys()


Out[10]:
['ombh2', 'omch2', 'w0', 'ns', 'ln10As', 'H0', 'Neff', 'z']

In [11]:
l = len(emu.scale_bin_centers)
idx = 0
params = {pname:pval for pname, pval in zip(emu._ordered_params.iterkeys(), emu.x[idx*l,:-1])}

In [12]:
for k,v in params.iteritems():
    print k,'\t'*5,v


Neff 					2.91875
H0 					63.36569
w0 					-0.8165972
omch2 					0.1140598
ln10As 					3.092918
ns 					0.975589
ombh2 					0.02268325

In [13]:
wp = emu.emulate_wrt_r(params, emu.scale_bin_centers)[0]
plt.plot(emu.scale_bin_centers, wp, label = 'Emu')
l = len(emu.scale_bin_centers)
plt.plot(emu.scale_bin_centers, emu.y[(idx)*l:(idx+1)*l]+emu.y_hat, label = 'Training')
plt.ylabel(r'$w_p(r_p)$')
plt.xlabel(r'$r_p \mathrm{[Mpc]}$')
plt.loglog();
plt.xscale('log')
plt.legend(loc='best')


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-13-fb0ced4f21ba> in <module>()
----> 1 wp = emu.emulate_wrt_r(params, emu.scale_bin_centers)[0]
      2 plt.plot(emu.scale_bin_centers, wp, label = 'Emu')
      3 l = len(emu.scale_bin_centers)
      4 plt.plot(emu.scale_bin_centers, emu.y[(idx)*l:(idx+1)*l]+emu.y_hat, label = 'Training')
      5 plt.ylabel(r'$w_p(r_p)$')

/u/ki/swmclau2/.conda/envs/hodemulator/lib/python2.7/site-packages/pearce/emulator/emu.pyc in emulate_wrt_r(self, em_params, r_bin_centers, gp_errs)
    669             del ep['z']
    670         elif 'z' not in self.fixed_params:
--> 671             raise ValueError("Please specify z in emulate_wrt_z")
    672         out = self.emulate_wrt_r_z(ep, r_bin_centers, z_bin_centers, gp_errs)
    673 

ValueError: Please specify z in emulate_wrt_z

In [ ]:
print wt
binlen = 19 for b in xrange(5): plt.plot(emu.scale_bin_centers, emu.y[b*binlen:(b+1)*binlen]+1) params = dict(zip(emu._ordered_params.keys(), emu.x[b*binlen+1,:-1])) pred = emu.emulate_wrt_r(params, emu.scale_bin_centers)[0] plt.plot(emu.scale_bin_centers, 1+pred, ls = '--') #plt.ylim([0.9, 1.5]) plt.xscale('log') plt.legend(loc='best') plt.xlim([0.1, 4]) plt.title(r'$w(\theta)$ Training vs. emulator') plt.xlabel(r'$\theta$') plt.ylabel(r'$w(\theta)$') plt.show()

In [ ]:
from pearce.mocks import compute_prim_haloprop_bins, cat_dict

In [ ]:
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[a]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!

cat.load_catalog(a, particles = True)
#halo_masses = cat.halocat.halo_table['halo_mvir']

In [ ]:
cat.load_model(a, 'hsabRedMagic')

In [ ]:
binno = 1
params = {pname:val for pname, val in zip(emu._ordered_params.iterkeys(), emu.x[binno*binlen,:-1])}
cat.populate(params)

wt = cat.calc_wt(theta_bins, do_jackknife=False,n_cores=1)

plt.plot(emu.scale_bin_centers,1+wt) plt.plot(emu.scale_bin_centers, 1+10**emu.y[binno*emu.scale_bin_centers.shape[0]:(binno+1)*emu.scale_bin_centers.shape[0]]) #plt.yscale('log') plt.loglog() plt.xlim([2.25/60, 275/60]) plt.ylim([0.8, 4.0])

In [ ]:
theta_bins = np.logspace(np.log10(2.5), np.log10(250), 20)/60
tpoints = (theta_bins[1:]+theta_bins[:-1])/2

In [ ]:
fig = plt.figure(figsize=(45,14))


emulation_point = [('f_c', 0.233), ('logM0', 12.0), ('sigma_logM', 0.333),
                    ('alpha', 1.053),('logM1', 13.5), ('logMmin', 12.033)]

em_params = dict(emulation_point)

em_params.update(fixed_params)
del em_params['z']

fixed_params2 = {'mean_occupation_satellites_assembias_param1':0.0,
                'mean_occupation_centrals_assembias_param1':0.0,
                'disp_func_slope_satellites':1.0,
                'disp_func_slope_centrals':1.0}

for idx, (param, bounds) in enumerate(emu._ordered_params.iteritems()):
    if param == 'r':
        continue
    wt_vals = []
    new_em_params = {}
    new_em_params.update(em_params)
    new_em_params.update(fixed_params2)
    for v in np.linspace(bounds[0], bounds[1], 6):
        new_em_params[param] = v
        wt_vals.append(emu.emulate_wrt_r(new_em_params, tpoints))
    wt_vals = np.array(wt_vals)
    
    pal = sns.cubehelix_palette(wt_vals.shape[0], start=idx, rot=0.3,\
                            dark=0.0, light=.60,reverse = True)
    #sns.palplot(pal)

    sns.set_palette(pal)

    #sns.set_style("darkgrid", {"axes.facecolor": "0.85"})
    plt.subplot(5,2,idx+1)

    for color, wt, v in zip(pal, wt_vals,np.linspace(bounds[0], bounds[1], 6) ):
        plt.plot(tpoints, 1+wt[0,:], color = color, label = r'%s = %.1f'%(param,v) )
    #plt.loglog()
    plt.xscale('log')
    plt.legend(loc='best')
    plt.xlim([0.1, 4])
    plt.title(r'$w(\theta)$ variance by %s'%param)
    plt.xlabel(r'$\theta$')
    plt.ylabel(r'$w(\theta)$')
plt.show()

In [ ]: